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Filament based intracellular transport involves the collective action of molecular motor proteins. 
Experimental evidences suggest that microtubule (MT) filament bound motor proteins such as 
kinesins weakly interact among themselves during transport and with the surrounding cellular en¬ 
vironment. Motivated by these observations we study a driven lattice gas model for collective 
unidirectional transport of molecular motors on open filament, which incorporates the short-range 
interactions between the motors on filaments and couples the transport process on filament with 
surrounding cellular environment through adsorption-desorption Langmuir (LK) kinetics of the mo¬ 
tors. We analyse this model within the framework of a Mean Field (MF) theory in the limit of weak 
interactions between the motors. We point to the mapping of this model with the non-conserved 
version of Katz-Lebowitz-Spohn (KLS) model. The system exhibits rich phase behavior with va¬ 
riety of inhomogeneous phases including localized shocks in the bulk of the filament. We obtain 
the steady state density and current profiles and analyse their variation as function of the strength 
of interaction. We compare these MF results with Monte Carlo simulations and find that the MF 
analysis shows reasonably good agreement as long as the motors are weakly interacting. We also 
construct the non-equilibrium MF phase diagram. 


Motor proteins play a crucial role in various intracellular processes. Their collective action is involved in host of 
functions ranging from intracellular trafficking to the organization of the complex cytoskeletal network [I]. Individual 
motor proteins such as kinesin and myosin hydrolyze ATP to actively move along the cellular filaments [TJ [2]. While 
single molecule experiments have shed light on the mechano-chemical properties and functioning of individual motor 
proteins [2|, it has been established that many of the cellular functions such as cell division and intracellular transport 
is achieved by the collective effect of many motors. Many motors working in unison not only generates large forces 
that are required for many cellular processes but it is also responsible for long distance regulated transport inside 
the cell 0. Thus it is of considerable interest to understand how the interaction between individual motor proteins 
affects the collective transport within the cell. 

Experimental studies on transport of kinesin motors on MT filament suggests that the individual motors interact 
weakly with their neigbours on the MT and that this attractive interaction manifests as a propensity of the motors 
to cluster during their transport on cellular filaments mg. It has also been widely observed that spatial organization 
of the motors and the cellular cargo that they transport on the filaments is sensitive to the surrounding cellular 
environment. For example, experiments done to study the spatial distribution of pigment granules in extracts of 
melanophore cell has illustrated that altering the cellular environment by specific biochemical means affects the 
binding affinity of the motors to the filaments and which in turn can lead to spatial reorganization of the pigment 
granules 000 . Further since cellular cargoes carried by the motors are loaded and offloaded at the filament ends, the 
role of the boundary processes at the filament ends also play a significant role in determining the spatial organization 
of the motors on individual filaments 0. 

In light of these experimental finding, we study a minimal lattice gas model for unidirectional transport on a single 
filament which takes into account the weak interaction between the motors and their interaction with the surrounding 
environment at the microscopic scales within a simplified framework. The proposed framework provides a natural 
means to qualitatively understand some of the aspects of collective behaviour of motors and their spatial organization, 
that emerges at macroscopic cellular scales. We study this model within a Mean Field (MF) theory and compare the 
results with Monte Carlo simulation results. 

The minimal model that we analyse, is a variant of the totally asymmetric exclusion process (TASEP) and belongs 
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FIG. 1: (color online) Schematic diagram of the process. 


to the broad class of one-dimensional drive diffusive system [9.. Driven diffusive systems in one dimension have been 
studied extensively from a theoretical perspective [MS] and they exhibit a variety of interesting features and rich 
phase behaviour which includes boundary induced phase transitions jT2til6] . Such boundary induced phase transitions 
in these driven systems are in contrast to one dimensional equilibrium systems where the boundaries do not affect 
the macroscopic phase of the system [12]. These models have proved to be particularly useful in modeling variety 
of processes such as transport across biomembrane channels HU, cellular cargo transport nans mim stochastic 
pumps pT} 122 ] . fungal growth [221 ED as well as the collective dynamics of interacting molecular motors 1251127] . 

The cytoskeletal filament on which the molecular motors move is represented as a one-dimensional lattice of length 
L. The motors bound to the filament are the ’particles’ in this description. The total number of lattice sites is N so 
that the lattice spacing e = L/N. The occupation number of a lattice site i is denoted by rii which assumes values 
rii = 1 if occupied and rq = 0 otherwise. Particles hop unidirectionally to the right on this lattice with excluded 
volume interactions preventing two particles from occupying a given site as in the usual TASEP model as long the 
particles are free (not part of any cluster). The short range attractive interaction between the motors is incorporated 
by including next-nearest neighbour (NNN) interaction. Due to the attractive interaction, the motors exhibit an 
enhanced tendency to form a cluster while the propensity to detach from a cluster is diminished. We include this 
effect by choosing a hop rate for a free particle (not part of a cluster) to jump into a cluster as q , with q > 1 while 
its hop rate out of motor cluster is chosen to be r = ^ [27]. Otherwise if a free particle hops to right and remains free 
then the hop rate is 1. We also choose a hop rate of 1 if a particle which is part of cluster jumps into another cluster. 
At the filament ends, particles enter the left boundary site with a rate a if the particle does not join a cluster, i.e; if 
the neighbouring site on the right of the boundary site is unoccupied, while the entry rate is qa if the particle enters 
the boundary and the adjacent site is occupied by a particle. For the boundary site at the filament end, the particles 
leave the lattice with a rate /3 as long the neighbouring site to the left of the right boundary site is unoccupied, 
while the the exit rate of particles is /3/q if the particle at boundary site is part of a cluster. The effect of embedding 
environment is incorporated by allowing for stochastic attachment and detachment process of motors on the filaments. 
This process breaks the particle conservation of the motors in the bulk of the filament. The introduction of Langmuir 
kinetic process of attachment and detachment is done such that a particle in the bulk can detach from a lattice site 
with a rate ujd while a particle could attach to a vacant lattice site with rate uja- We restrict ourselves to the regime 
where the boundary processes of particle input and output compete with the bulk process of stochastic attachment 
and detachment. Thus we define bulk attachment rate Qa = uaN and bulk detachment rate Qd = wdN and hold 
the bulk attachment and detachment rates fixed as N oc HSUI6]. 

The various processes are summarized in Fig. 1. We normalize the the total length of the lattice L to 1 . The 
continuum Mean Field (MF) evolution equation of mean density of particles is obtained by taking the expectation 
value of site occupation numbers, invoking the usual procedure of factorizing the two point correlators arising out of 
different combinations of site occupation numbers and then expressing the resultant equations in terms of continuum 
variable of relative position x along the filament lattice [muni eh]. The continuum MF evolution equation in the 
bulk reads as, 
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where the current has the form 


J{p) = p{ i - p)[(i - pf + p 2 ] + (q + ^)p 2 (! - p) 2 - ( 2 ) 

Here we have displayed terms unto the leading order in 5. 

The corresponding equation for steady state density profile p(x) is obtained by simply setting the time derivatives of 
p in Eq. 0to zero. For satisfying the left boundary condition at x = 0, the density profile has to satisfy the conditions 
p(0) = a and for density profile satisfying the right boundary condition, p(l) = 1 — /3. 

It is worthwhile pointing out all the different limits of this model. In the absence of next nearest neighbour (NNN) 
interaction and Langmuir kinetics (LK) of particle attachment and detachment, the dynamics of the model reduces 
to the usual TASEP. In this limit, which is obtained by setting q = 1 and 1a = 0, we obtain steady state MF 
equation which corresponds to MF steady state solution for TASEP. For the situation when LK is present, while 
the NNN interaction is absent, the steady state equation reduces to steady state equation for the model discussed in 
ns esi, wherein the steady state allows for phase coexistence and shocks localized in the bulk. Finally when the the 
Langmuir kinetics is absent while NNN interaction is included, then the MF steady state solution is same as the steady 
state density solution discussed in [28] . In fact in this limit, the model maps on to Katz- Lebowitz- Spohn (KLS) 
model [29] [30]. For KLS model particle hopping is dependent on NNN and is defined by the rules: 1011 —»• 1101, 
1010 —)> 1100, 0010 —)> 0100 and 0011 0101, where 1 corresponds to a particle occupancy at a lattice site while 0 

corresponds to a vacancy, and where these processes occur in general with different rates nans. In the absence of 
LK, the model discussed in this letter exactly maps to the dynamics of the KLS model, by replacing particles of the 
KLS model by vacancies and vacancies of KLS model by particles for our case. While for KLS model it has been 
shown that MF analysis breaks down for strong NNN interaction, we restrict ourselves to the situation when the NNN 
interaction is weak , as is the case for motor interactions on cytoskeletal filaments. In this limit the MF picture can 
fairly accurately describe the steady state behaviour as confirmed by Monte Carlo simulations. Further, unlike KLS 
model, here the interplay of the bulk translational dynamics on the lattice with the Langmuir Kinetic processes gives 
rise to qualitatively new scenario, wherein we find that phase coexistence and localized shocks in the bulk which are 
also controlled by attachment-detachment rates associated with LK process apart from the NNN coupling strength 
and the boundary entry and exit rates. Thus the model discussed in this letter can be considered as a non-conserved 
version of the KLS model. 

In order to simplify the analysis and construct the entire MF phase diagram in the weak motor interaction limit, 
we restrict our analysis for a special choice of Langmuir kinetic processes wherein the attachment rate tip is equal to 
the detachment rate Qa- For this special case, when Qd = £Ia = the steady state differential equation for density 
is, 


(1 - 2 p) 
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Like in case of TASEP, one of the possible solution of Eq.(|3| is the homogeneous solution, p{pc) = This corresponds 
to the maxima of current in the lattice. The corresponding expression for the current in this maximal current (MC) 
phase is Jm = § + p^(g + l/g) an d the density profile becomes independent of the boundary densities. Unlike TASEP, 
the maximal current in the lattice for this case is Jm < 1/4 for q ^ 1. Apart from this homogeneous solution for the 
density, the other possible solution is obtained by integrating Eq. and is in form of a form a cubic equation in p, 

P 3 - 2 P 2 ~ + ~]^ x + C) = 0 ( 4 ) 


where A = 2(q + 1 /q) — 4 and C is a constant which is fixed by appropriate boundary condition. Using the mean-field 
boundary conditions in the steady state, p(0) = a and p(l) = 1 — /3, we can determine C from Eq. ©• Substituting 
C back in Eq. [4] we can then solve for the cubic equation to determine the physically relevant solutions which will 
give the density profile p(x). This is then used in Eq. [2j to calculate the current in the different phases. When the 
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FIG. 2: (color online) Average density p(x) and current j(x) profiles for parameters a = 0.2, ^ = 0.1,0 = 0.3 and (a) q = 1.0 
(b) q — 1.2 and (c) q— 1.6. Profiles calculated analytically and from Monte Carlo simulations. 


density profile satisfies the boundary condition at the left boundary at x = 0 , p(0) = a is substituted in Eq. Q 
to obtain the corresponding expression for C and the entire density profile in the range 0 < x < 1 is subsequently 
obtained by solving Eq. This corresponds to the Low Density (LD) solution. Similarly the High Density (HD) 
solution for the density profile is obtained by by using the boundary condition p(l) = 1 — /?. It is important to note 
that Eq. [3] is a first order differential equation which has to satisfy two boundary conditions. Therefore, the problem 
is overdetermined and the two solutions that respects the two boundary conditions will in general not match and give 
rise to a shock in the bulk. The position of the shock is determined by the continuity of the current J for the two 
different solutions. In general three different density profiles may coexist on the filament lattice, e.g; the LD, HD and 
the MC phase. For any two phases to coexist, the current J corresponding to the two different phases at a particular 
location x = x Q on the lattice must be equal for the domain range 0 < x < 1, x Q being the position of the domain 
wall separating the two phases. While the density profile is continuous across a domain wall separating a LD or HD 
phase region with a MC phase region, the density change across a domain wall separating a LD phase with HD phase 
is discontinuous and results in a shock localized in the bulk (Fig. [2]). The location of the domain wall coincides with 
the position where the current is maximum, the current increasing and decreasing monotonically with x on either 
side of it. With increasing interaction (increasing g), the domain wall shifts to the left. For very strong interaction, 
the domain wall disappears as the system enters into the HD phase, the density profile being larger than 1/2. As 
expected, the mean field results no longer hold for large q. 

In general the possible phases in the bulk of the filament are the pure LD, HD and MC phases, two-phase coexistence 
of LD with MC (LD-MC), HD with MC (HD-MC) and LD with HD (LD-HD), and three-phase coexistence of LD 
with MC and HD ( LMH). LD-HD phase is characterized by density discontinuity and shock localized in the bulk of 
the lattice. The phase boundaries separating the different possible phases are determined by setting the domain wall 
position to left or the right end of the filament lattice. We now discuss the conditions which determines the MF phase 
boundaries that separates the different possible phases. 

Phase coexistence boundary between LD and LD-HD : Starting from Eq. the density profile corresponding to LD 
phase is obtained by using the boundary condition p(0) = a. The corresponding current profile J{x) is obtained by 
using Eq. <©• Similarly the current solution for HD phase is obtained by using the boundary condition p(l) = 1 — /?. 
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FIG. 3: (color online) Mean field phase diagram: (a) q = 1.2, Q = 0.3, (b) q = 2.0, Q = 0.3. 


To obtain the phase boundary, the current solution corresponding to LD phase Jld is equated with the current 
solution corresponding to HD phase Jhd at x = 1, i.e; Jld{% = 1) = Jhd{x = 1). 

Phase coexistence boundary between HD and LD-HD : For obtaining the phase boundary, the current solution 
corresponding to LD phase Jld is equated with the current solution corresponding to HD phase Jhd at x = 0, i.e; 
Jld(x = 0) = Jhd{x = 0). 

Phase coexistence boundary between LD and LD-MC : For obtaining the phase boundary in this case, the density 
solution for the low density phase pld is equated to 1/2 at the right boundary, i.e; Pld(% = 1) = The resultant 
phase boundary has the equation of a straight line. 

Phase coexistence boundary between MC and LD-MC : Here the phase boundary is determined by equating the 
density solution for the LD phase phd to 1/2 at the left boundary of the filament, i.e; Pld{% = 0) = This is simply 
an equation of straight line a = 1/2. 

Phase coexistence boundary between MC and HD-MC: For obtaining the phase boundary in this case, the density 
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FIG. 4: (color online) Steady state current as a function of q for a = 0A,/3 = 0.4, Q = 0.0. Points are simulations with error 
bars; line shows fit to ~ 1 /q for large q. 


solution for the high density phase phd is equated to 1/2 at the right boundary, i.e; Phd{% = 1) = i- Thus the 
resultant phase boundary is a straight line (3 = 1/2. 

Phase coexistence boundary between HD and HD-MC : For obtaining the phase boundary in this case, the density 
solution for the high density phase phd is equated to 1/2 at the left boundary, i.e; Phd{% = 0) = The resultant 
phase boundary is a straight line. 

Phase coexistence boundary between LD-M-HD and LD-HD : This is obtained by requiring that the density solution 
for the high density phase Phd, low density phase pld are simultaneously equated to 1/2 at a location x in the bulk 
of the lattice. 

In Fig. 3 we display the MF phase diagram for different values of q. When q is sufficiently small, the topology of 
the phase diagram is very similar to that of TASEP-LK phase diagram as in Ref. [16, . With increasing value of q the 
location of the phase boundaries changes. Further unlike TASEP-LK, the phase boundary curve separating LD-HD 
phase with HD and the phase boundary curve separating LD-HD phase with LD is not a straight line. For high values 
of q corresponding to strong NNN interaction, the MF predications significantly deviates from the actual density and 
current profiles and the corresponding phase boundaries. In the limit of q » 1, the MF picture completely breaks 
down. In the absence of LK process, the current in the lattice in this limit is controlled solely by the slowest process in 
the lattice, i.e; the rate at which the particles leave the cluster. In the absence of Langmuir kinetics, the corresponding 
steady state current in the lattice J ~ 1/g, for large q. Fig.4 shows this variation. 

In summary we have studied a driven lattice gas model for the transport of molecular motors which not only interact 
with each other but also with cellular environment through adsorption-desorption Langmuir kinetics. We showed that 
the resulting phase diagram is extremely rich, with a variety of inhomogeneous phases, including localized shocks in 
the bulk. Our mean field results match reasonably well with Monte Carlo simulations for weak interactions. The 
model maps to the non-conserved version of KLS model. For strong interactions between particles, the simple MF 
theory breaks down. 
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